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Abstract 

The paper deals with polynomial interpolation, least-square approximation and cu- 
bature of functions defined on the rectangular cylinder, K = D x [—1,1], with D the 
unit disk. The nodes used for these processes are the Approximate Fekete Points (AFP) 
and the Discrete Leja Points (DLP) extracted from suitable Weakly Admissible Meshes 
(WAMs) of the cylinder. From the analysis of the growth of the Lebesgue constants, 
approximation and cubature errors, we show that the AFP and the DLP extracted from 
WAM are good points for polynomial approximation and numerical integration of func- 
tions defined on the cylinder. 

1 Introduction 

Locating good points for multivariate polynomial approximation, in particular interpolation, 
is an open challenging problem, even in standard domains. 

One set of points that is always good, in theory, is the so-called Fekete points. They are 
denned to be those points that maximize the (absolute value of the) Vandermonde determinant 
on the given compact set. However, these are known analytically only in a few instances 
(the interval and the complex circle for univariate interpolation, the cube for tensor product 
interpolation), and are very difficult to compute, requiring an expensive and numerically 
challenging nonlinear multivariate optimization. 

Admissible Meshes (shortly AM), introduced by Calvi and Levenberg in [8], are sets of 
points in a given compact domain K C M. d which are nearly optimal for least-squares approx- 
imation, and contain interpolation points that distribute asymptotically as Fekete points of 
the domain. This theory has given new insight to the (partial) solution of the problem of 
extracting good interpolation point in dimension d > 1. In all practical applications, instead 
of AM, people look for low-cardinality admissible meshes, called Weakly Admissible Meshes 
or WAM (cf. the recent survey [5]). 

The extremal sets of our interest are the Approximate Fekete Points (AFP) and Discrete 
Leja Points (DLP). As described in [U Ho] , AFP and DLP can be easily computed by using 
basic tools of numerical linear algebra. In practice, (Weakly) Admissible Meshes and Discrete 
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Extremal Sets allow us to replace a continuous compact set by a discrete version, that is "just 
as good" for all practical purposes. 

In this paper, we focus on AFP and DLP extracted from Weakly Admissible Meshes of the 
rectangular cylinder K = D x [—1,1], with D the unit disk. These points are then used for 
computing interpolants, least-square approximants and cubatures of functions defined on the 
cylinder. We essentially provide AFP and DLP of the cylinder that, to our knowledge, have 
never been investigated so far. 

2 Weakly Admissible Meshes: definitions, properties 
and construction 

Given a polynomial determining compact set K C M d or K C C d (i.e., polynomials vanishing 
there are identically zero), a Weakly Admissible Mesh (WAM) is defined in [8] to be a sequence 
of discrete subsets A n C K such that 

\\p\\K<C{A n )\\p\\ An , VpeF d n (K) (2.1) 

P^(fT) being the set of ci-variate polynomials of degree at most n on K, where both card(A n ) > 
iV := dim(P^(fr)) and C(A n ) grow at most polynomially with n, i.e. card(y4 n ) < c n s , for 
some fixed s G N depending only on K. When C(A n ) is bounded we speak of an Admissible 
Mesh (AM). We use the notation WfWx = sw P x eK \f(x)\ for / a bounded function on the 
compact K. 

WAMs enjoy the following ten properties (already enumerated in |3] and proved in [8]): 
PI: C(A n ) is invariant under affine mapping 

P2: any sequence of unisolvent interpolation sets whose Lebesgue constant grows at most 
polynomially with n is a WAM, C(A n ) being the Lebesgue constant itself 
P3: any sequence of supersets of a WAM whose cardinalities grow polynomially with n is a 
WAM with the same constant C(A n ) 

P4: a finite union of WAMs is a WAM for the corresponding union of compacts, C(A n ) being 
the maximum of the corresponding constants 

P5: a finite cartesian product of WAMs is a WAM for the corresponding product of compacts, 
C(A n ) being the product of the corresponding constants 

P6: in C d a WAM of the boundary OK is a WAM of K (by the maximum principle) 

P7: given a polynomial mapping 7r s of degree s, then 7i s (A ns ) is a WAM for ir s (K) with 

constants C(A ns ) (cf. [3j Prop. 2]) 

P8: any K satisfying a Markov polynomial inequality like ||Vp||jc < Mn r \\p\\k has an AM 
with 0{n rd ) points (cf. [SJ Thm.5]) 

P9: least-squares polynomial approximation of / G C(K): the least-squares polynomial L A „f 
on a WAM is such that 

||/ - L A J\\ K < C(A n Wcavd(A n ) min{||/ - p\\ K , p G F d (K)} 

(cf. Thm.l]) 

P10: Fekete points: the Lebesgue constant of Fekete points extracted from a WAM can be 
bounded like A n < NC(A n ) (that is the elementary classical bound of the continuum Fekete 
points times a factor C(A n )); moreover, their asymptotic distribution is the same of the 
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continuum Fekete points, in the sense that the corresponding discrete probability measures 
converge weak-* to the pluripotential equilibrium measure of K (cf. [3j Thm.l]). Pluripotential 
theory has been widely studied by M. Klimek in the monograph [TU], to which interested 
readers should refer for more details. 

It is worth noticing that in the very recent papers [Tlj 113] , the authors have provided new 
techniques for finding admissible meshes with low cardinality, by means of analytical trans- 
formations of domains. 



Examples of WAMs can be found in jH [5] . Here, we simply recall some one dimensional and 
two dimensional WAMs. 

1. The set 

C n = {cos(kir/n), k = 0,...,n} 

of n + l Chebyshev-Lobatto points for the interval I = [—1, 1], is a one-dimensional WAM 
of degree n with C(A n ) = 0(log n) and card(C n ) = n + 1. This follows from property 
P2. 

2. The set Pad„, n > of the Padua points of degree n of the square Q = [— 1, l] 2 is the 
set defined as follows (cf. [2]) 

Pad n = {x kd = (£ k ,r h ), 0<A;<n, 1 < j < Lf J + 1}, (2.2) 

where 

, ( cos ^tttt, k even 

kit n+1 
£ fc = cos — , Vj = \ ( 2 -3) 

I cos ^rrTC, k odd 

V n+l ' 

Notice that here we refer to the first family of Padua points. Pad n is then a two- 
dimensional WAM with C(Pad n ) = C(log 2 n) and card(Pad n ) = (n + l)(n + 2)/2. This 
is a consequence of property P2 since, as shown in [2], the Padua points are a unisolvent 
set for polynomial interpolation in the square with minimal order of growth of their 
Lebesgue constant, i.e. C(log 2 n). 

3. The sequence of polar symmetric grids A n = {(rj cos 9 j,ri sin 6 j)} with the radii and 
angles defined as follows 

{n,9j)i,j = {cos(ivr/ri), < i < n} x < j < nj (2.4) 

are WAMs for the closed unit disk D = {x : ||x|| 2 < 1}, with constant C(A n ) = 
(9(log 2 n) and cardinality card(A n ) = (n + l) 2 for odd n and card(A n ) = n 2 + n + 1 for 
even n (cf. [6, Prop. 1]). Moreover, since these WAMs contain the Chebyshev-Lobatto 
points of the vertical diameter 6 = n/2 only for n odd (whereas it always contains the 
Chebyshev-Lobatto points of the horizontal diameter 6 = 0), and thus is not invariant 
under rotations by an angle 7r/2. Hence in order to have a WAMs on the disk invariant 
by rotations of tt/2, we have to modify the choice of radii and angles in (2.4) as follows 

{J7T ] 
0<j<n + l>,n even (2.5) 

In this way the obtained WAM is now invariant with card(A n ) = (n + l) 2 also for n 
even. 
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2.1 Three dimensional WAMs of the cylinder 

We restrict ourselves to the rectangular cylinder with unitary radius and height the interval 
[-1,1], that is K = D x [—1, 1], where as above, D is the closed unit disk. 

We considered two meshes: the first one uses a symmetric polar grid in the disk D and 
Chebyshev-Lobatto points along [—1,1]; the second one uses Padua points on the (x, z) plane 
and equispaced points along the circumference of D. 

2.1.1 The first mesh: WAM1 

We consider the set 

A n = { (r , cos 9j , r { sin Oj ,z k )} 
with — 1 < Ti < 1, < % < n and < 9j < 7T, < j < n that is 

r/ f /«7r\^^.^l f ^t^, < j < n + 1, neven 

I VW J I ^>0<j<n, nodd 

x <( cos I — ),0</c<n 



The cardinality of A n , both for n even and n odd, is (n + l) 3 . Indeed, let us consider first 
the case of n even. The points on the disk, subtracting the repetitions of the center, which 
are n + 2 — 1, are (n + l) 2 . All these points are then multiplied by the corresponding n + 1 
Chebyshev-Lobatto points along the third axis z, giving the claimed cardinality. 
When n is odd, there are no coincident points, thus we have (2n + 2){n + l)/2 = {n + l) 2 
points on the disk. Then, considering the n + 1 Chebyshev-Lobatto points along the third 
axis, we get the claimed results. 

Finally, the set A n so defined, is a WAM since it is the cartesian product of a two dimensional 
WAM (the points on the disk) and the one dimensional WAM of the Chebyshev-Lobatto 
points. The property P5 gives the constant C(A n ) = (9(log 3 n) (see Figure [3] for the case 
n = 5). 

2.1.2 The second mesh: WAM2 

This discretization is obtained by taking the Padua points Pad n on the plane (r, z), rotated 
n times along z-axis by a constant angle 9 — 7r/(n + 1). In this way, along the bottom 
circumference of the cylinder, we obtain In + 2 equispaced points. This is due to the fact that 
the points with coordinates (—1,0) and (1,0) are Padua points. In details, the mesh is the set 

A n = { (n cos 9j , Ti sin 9j ,z k )} 

with — 1 < Ti < 1, < i < n and < 6j < tt, < j < n, that is 

{(ri,9j,z k )}ij,k = j cos ,0 < i < n\ x j^->0 < j < n + lj 

f / kn \ A; odd when i is even 

x < cos ,0<A;<n + l , , .. 

{ \n + l J k even when 1 is odd 

This mesh has cardinality O {^^j ■ hi fact, when n is even, the points are ( ra+1 K n + 2 ) ^ + 1) from 
which we have to subtract the repetitions + 1) n, corresponding to the f + 1 Padua points 
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with abscissa x = counted n times. Then, the points so generated are (n 2 + n + 1) ( w + 2 ) . On 
the contrary, when n is odd, there are no intersections and so the total number of points is 

(n+lHra+2) ^ gee pjg ure |3 f or cage /J, = 5,6). 

We now prove that this mesh is indeed a WAM. To this aim, consider a generic polynomial 
of degree at most n defined on the cylinder p(x,y,z) = p(r cosd, r sin 8, z) := q(8,r,z). For 
a fixed angle 8, q is a polynomial of degree at most n in r, z, while it is a trigonometric 
polynomial in 8 of degree at most n for fixed values of (r, z). Since on the generic rectangle 
(r, z) we have considered the set of Padua points of degree n which is a WAM, say A\, hence, 
we can write 

\q(0,r,z)\ <ci(n)\\q(8,;-)\\ Al , 

where c\(n) does not depend on 8. Let \q(8,r* , z*)\ be the maximum. Considering now the 
equispaced angles, 8 k = 2/c7r/(2n + 2), < k < 2n + 1 i.e. the 2n + 2 equispaced points in 
[0, 2n[, say A 2 , then 

\q(lr*,z*)\ <c 2 (n)\\q(;r*,z*)\\ A2 . 
Passing to the maximum also on the left side, we have 

max \p(x,y,z)\ = \\p\\ K < ci(n)c 2 (n)\\q\\ An 

that is \\p\\k < C( n )lblU„) where C(n) is indeed C(A n ) = 0(\og 3 n), showing that this 
discretization is a WAM for the cylinder K. 
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Figure 1: Left: Padua points for n = 5 on the (r, z)-plane. We used different markers for 
the two Chebsyhev meshes of Padua points. Right: Projections on the disk of WAM2 points 
(obtained by n + 1 rotations of Padua points around the z axis). Notice that the points 
indicated with the small triangles are on different levels. See also Figure [?]. 



3 Computation of AFP and DLP 

As discussed in [1] , the computation the AFP and DLP, can be done by a few basic linear alge- 
bra operations, corresponding to the LU factorization with row pivoting of the Vandermonde 
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Figure 2: Left: Padua points for n = 6 on the (r, z)-plane. We used different markers for 
the two Chebsyhev meshes of Padua points. Right: Projections on the disk of WAM2 points 
(obtained by n + 1 rotations of Padua points around the z axis). Notice that the points 
indicated with the small triangles are on a different levels. 



matrix for the DLP, and to the QR factorization with column pivoting of the transposed Van- 
dermonde matrix for the AFP (cf. [15]). For the sake of completeness, we recall these two 
Matlab-like scripts used in [151 H] f° r computing the AFP and DLP, respectively, 
algorithm AFP (Approximate Fekete Points): 

• W= (V(a,p)Y; b = (1,...,1)' G C^; w = W\b; ind = f ±nd(w ^ 0); £ = a(ind) 
algorithm DLP (Discrete Leja Points): 

• V = V(a,p); [L,U,a] = LU(V, "vector"); ind = cr(l,...,A); £ = a(ind) 

In Figures [4]-[5j we show the AFP and DLP extracted from the WAM1 and WAM2 for 
n = 5. 

In the above scripts, V(a,p) indicates the Vandermonde matrix at the WAM a using the 
polynomial basis p, that is the matrix whose elements are Pj(a,i), 1 < j ' < N, 1 < i < 
card(v4 n ). The extracted AFP and DLP are then stored in the vector 

Remark 1. In both algorithms, the selected points (as opposed to the continuum Fekete 
points) depend on the choice of the polynomial basis. But in the second algorithm, which is 
based on the notion of determinant (as described in [H §6.1]), the selected points also depend 
on the ordering of the basis. In the univariate case with the standard monomial basis, it is 
not difficult to recognize that the selected points are indeed the Leja points extracted from 
the mesh (cf. [H [H] and references therein). 

Remark 2. When the conditioning of the Vandermonde matrices is too high, and this happens 
when the polynomial basis is ill-conditioned, the algorithms can still be used provided that 
a preliminary iterated orthogonalization, that is a change to a discrete orthogonal basis, is 
performed (cf. [31 HI [15]). This procedure however only mitigates the effect of a bad choice of 
the polynomial basis. Consequently, whenever is possible, is desirable to use a well-conditioned 
polynomial basis. 

A suitable basis for the forementioned rectangular cylinder K, is the set of polynomials intro- 
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Figure 4: WAM1 of the cylinder for n = 5 and the corresponding extracted points. Left: 56 
Approximate Fekete Points. Right: 56 Discrete Leja Points. 
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Figure 5: WAM2 of the cylinder for n = 5 and the corresponding extracted points. Left: 56 
Approximate Fekete Points. Right: 56 Discrete Leja Points. 
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duced by J. Wade in [18] : 

C jjk>i (x,y,z) :=U k {0j,k;x,y)fi- k (z), i = 0, . . . ,n, k = 0, . . . ,i, j = 0, . . . , k (3.6) 

where 

• ®h k = k+i'i 

• Uk(9j t k] x,y) = U k (x cos(6 ) j i / c ) + ysm(9j ik )) is the Chebyshev polynomial of the second 
kind which is an orthonormal basis for the disk w.r.t. the measure u(x,y) = -; 

• Tj(z) is the j-th orthonormal Chebyshev polynomial of the first kind, i.e. Tj(z) = 
y/2Tj(z), w.r.t. the measure u{z) = (1 — z 2 )~^. 

As discussed in [15] . the basis Cj yk ,i{x,y, z) is orthonormal for the space of orthogonal 

polynomials on K w.r.t. the weight function ^(1 — z 2 ) 2 . This basis plays also an important 
role in the construction of the discretized Fourier orthogonal expansions on the disk D and 
the unitary rectangular cylinder. Moreover, this turns out to be well-conditioned. Indeed, as 
outlined in [TH] , if we consider the Radon projection of a function / on W 1 



Re(f;t)= [ f{x)dx, teR, 6e 
J<x,e>=t 



3d— 1 



with § d 1 the unit sphere in Mr. The discretized Fourier expansion F 2m of Re(f; t) for m > 0, 
(x,y) belonging to the unit disk and z G [—1, 1] 

2m 2m 2m— 1 

F 2m (f)(x,y,z) = ^(/(■'■' z fc,2m;cos(^- 2m )))7^ ijifc (x,|/, z) 

1=0 j=l k=0 



where cp u 



2un 
2m+l> 

2m n 



Ti,j, k {x, y, z) = 4m 2^ 2m + l) z2( l + ^ sin( ^ + l )9j,2m)Ui{cos{ai{x, y)))T n „i(z k ,2m)T n -i(z) 



n=0 1=0 



and (Ti(x, y) = arccos(x cos <fii+y sin0j), has sup norm that grows as ||F2 m ||oo ~ m(log(m+l) 2 ), 
i.e. nearly the optimal growth. 



4 Approximation and cubature on the cylinder 
4.1 Interpolation and least-square approximation 

The interpolation polynomial g n ( x ) of degree n of a real continuous function / defined on the 
compact K C M 3 , can be written in Lagrange form as 

N 

^^Mxj.xeif (4.7) 
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where N = dim(P„), a 3 are the AFP or the DLP extracted from the WAMs and lj indicates 
the jth elementary Lagrange polynomial of degree n. Let 1 = (Zi(x), . . . , Zjv( x )) be the (row) 



vector of all the elementary Lagrange polynomials at a point x, p the vector of the basis (3.6) 
and a = (ai, . . . , a^v), then we can compute 1 by solving the linear system 

1* = Wp\ W = (V^p)- 1 )*. 

The interpolation operator C n : C(K) — > P^, with C(K) equipped with the sup norm, that 
maps every / E C{K) into the corresponding polynomial q = Y^iLi ^(')/( x «) m Lagrange 
form, is a projection having norm 



N 

\C n \\oo = maxV]|Z i (x)| := A n (4.8) 

i=l 



where A n is the well-known Lebesgue constant. When the interpolation points are the true 
Fekete points, the Lebesgue constant satisfies the upper bound 

A n = max ||Wp*||i < N 

xeif 

since ||Zj|| < 1. 

Thanks to property P10 of WAMs we can say more. Indeed, when the Fekete points are 
extracted from a WAM, < C(y4 n )||/ J || yln Vj (cf. [3 §4.4]), from which the following 

upper bound holds 

K<NC{A n ), 

with C(A n ), the same constant in definition of a WAM, which depends on A n . In the numerical 
experiments that we will present in the next section, we will observe that the above upper 
bound is a quite pessimistic overestimate. 

Another natural application of such a construction, is the least squares approximation of a 
function / e C{K). Given a WAM A n = {a 1; . . . , a M }, M > N, with N = dim(p3) and 
P — {Pi; • • • iVn} a basis for P^, let us consider the orthonormal basis q w.r.t. the discrete 
inner product (/, g) = Yli=i /( a «)5'( a i) which can be obtained from the basis p by multiplying 
by a certain transformation matrix P, i.e. q=Pp. 

The least squares operator of / at the points of a WAM A n can then be written as 

N I M \ M 

j=\ \i=i j i=i 

where gfj(x) = J2j=i Qji^Qji 3 ^)- Letting g = (g u . . .,g M ) f , it follows that g = QP*p, where 
the matrix Q is a numerically orthogonal (unitary) matrix, i.e. Q l Q = I, and Q = V(a, q) = 
V(a, p)P. Notice that, the transformation matrix P and the matrix Q are computed once 
and for all for a fixed mesh. 

The norm of the operator, that is its Lebesgue constant, is then 

M 



V|^(x)| = max||gP* P (x) 



max 

8=1 
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In [8], it is observed that three-dimensional WAMs so defined can be used as discretization of 
compact sets K C M 3 for the computation of good points for least-square approximation by 
polynomials. Indeed, using property P9, in [HI Th. 2] the authors proved the following error 
estimates for least-squares approximation on WAMs of a function / e C(K) 

||/ -L An {f)\\ K < (l + C(A0(1 + VcardtAJ)) min{||/ - p\\ K : p E P 3 } (4.9) 

where again C(A n ) depends on the WAM. 

This estimates says that if we could control the factor C(A n )(l + a/ card(v4„)), and this is 
possible when we use WAMs, then the approximation L^ n (f) e P 3 is nearly optimal. 

4.2 Cubature 

For a given function / : K C IR 3 — > K we want to compute 

1(f) = [ /(x)dx 

where dx is the usual Lebesgue measure of the compact set K. An interpolating cubature 
formula Cjv(/) that approximates /(/) can be expressed as 

N 
i=l 

where, in our case, the nodes Xj are the AFP or the DLP for K. Once we know the cubature 
weights Wi, the Cjv(/) gives an approximation of 1(f). The cubature weights can be determined 
by solving the moment system, that is 

N 

WjPiixj) = / Pi(x)dx, i = 1, . . . , N 

3=1 JK 

where Pi is the z-th element of the polynomial basis p. Hence, if V — (pi(xj)) and &j = 
J K Pi(x)dx, the nodes Xi and the weights Wi are provided by the AFP or DLP algorithm. 

5 Numerical results 

In this section we present the numerical experiments that we made for showing the quality 
of AFP and DLP on interpolation, approximation by least-squares and cubature of 6 test 
functions defined on the rectangular cylinder K = D x [—1,1]. The results are obtained by 
using both the AFP and the DLP extracted from both WAM1 and WAM2. 

In Table [T] we collect, for degrees n = 5, 10, 15,20,25,30, the values of the Lebesgue con- 
stant A n , those of the condition number (in the sup norm) K ra]1 for the associated Vandermonde 
matrix using the Wade basis for the WAM1 for the AFP. Due to hardware restrictions, the 
Lebesgue constant has been evaluated on a control mesh M n = A m with m = An for n < 20, 
m = 2n for n > 20. 

In Table [2] we present A n , K^a^ f° r the WAM2 on the AFP. In this case the control mesh 
M n = A m with m = 4n for n < 20, m = 3n for n < 25 and m = 2n for n > 25. 
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In Tables [3] and [4] we present the corresponding values for the DLP. 

In Table [5] we display the the norm of the least-square operator ||-Zm„|| on the WAM1 and 
WAM2, respectively. 

Note that the norm of the least-square operator ||ImJ| (norm computed at the points of 
the WAM A n ) has been computed after two steps of orthonormalization of the polynomial 
basis. 

The numerical results show that both A n and the condition number are smaller for 
WAM2 while, by oppositeon the contrary, the sup norm of the least-square operator turns out 
to be bigger for WAM2 than WAM1. Actually, for WAM2, ||LaJ| has a growth factor close 
to 2. 



n 


5 


10 


15 


20 


25 


30 


A„ 


17 


83 


208 


384 


849 


988 




18.2 


177 


384 


746 


1410 


2650 



Table 1: The Lebesgue constant and the condition number of the Vandermonde matrix con- 
structed using the Wade basis on the AFP for the WAM1. 



n 


5 


10 


15 


20 


25 


30 


A„ 


19 


76 


213 


427 


879 


1034 


K oo,2 


19.4 


115 


440 


705 


1540 


2380 



Table 2: The Lebesgue constant and the condition number of the Vandermonde matrix con- 
structed using the Wade basis on the AFP for the WAM2. 



n 


5 


10 


15 


20 


25 


30 


A„ 


30 


115 


350 


617 


1388 


2597 


K oo,l 


35.2 


247 


772 


1190 


3090 


5320 



Table 3: The Lebesgue constant and the condition number of the Vandermonde matrix con- 
structed using the Wade basis on the DLP for the WAM1. 
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n 


5 


10 


15 


20 


25 


30 


An 


30 


129 


349 


648 


1520 


2143 


Koo,2 


29.9 


176 


638 


782 


2310 


3940 



Table 4: The Lebesgue constant and the condition number of the Vandermonde matrix con- 
structed using the Wade basis on the DLP for the WAM2. 



n 


5 


10 


15 


20 


25 


30 


II^JI 


4.8 


10.2 


10.7 


21.1 


15.8 


22.6 


\\L An J 


7.2 


15.3 


32.8 


43.4 


85.9 


96.6 



Table 5: The sup-norm of the least-squares operator on WAM1 and WAM2, respectively. 



The interpolation and the cubature relative errors on the AFP and the DLP have been 
computed for the following six test functions: 



f2(x,y,z 

h(x,y,z 
fe(x,y,z 



0.75e- 



(!)x+l)~ 9y+l <3z+l 



10 10 



+ 0.75e 49 

+ 0.5e~ J ' '"4 ' _ .2e-( 9a: - 4)2 - ( '" / - 7,i - , ' i "-- : ' 1 

= ^(x - OA) 2 + (y- 0.4) 2 + (z- OA) 2 ; 
= cos(4(x + y + z)); 



1 + 16(x 2 + y 2 + z 2 ) 
^/{x 2 + y 2 + z 2 f; 
cos(x 2 + y 2 + z 2 ) . 



The function fx is the three-dimensional equivalent of the well-known Franke test function. 
The function /2 has a singular point into the cylinder K = D x [—1,1]. The functions fs and 
fe are infinitely differentiable. The function f± is the Runge function. The function f§ is a C 2 
function with third derivatives singular at the origin. 

All numerical experiments have been done on a cluster HP with 14 nodes. We used one of 
the nodes equipped with 2 processors quad core with 64Gb of RAM. 

In Figures [6] and [7] we display the interpolation, cubature and least-square relative errors on 
the AFP and DLP, up to degree n = 30, for the WAM1 and WAM2, respectively. The results 
shows that the AFP give, in general, smaller errors. Only the cubature errors on WAM2 are 
smaller for DLP than AFP. One reason is related to the values of the Lebesgue constants and 
the conditioning of the Vandermonde matrices that are smaller for AFP than DLP. 

Since WAM2 has a lower cardinality and that the results are more or less the same, such a 
mesh is more convenient from the point of view of efficiency and approximation order. As true 
values of the functions fa, i = 1, . . . , 6, we considered the value of fi on the control meshes 
used for computing the Lebesgue constants. As exact values of the integrals, we considered 
the values computed by the Matlab built-in function triplequad with the chosen tolerance 
depending on the smoothness of the function. For the smoothest functions fs and /g we used 
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the tolerance l.e — 12 while for the others l.e — 10. This choice allowed to avoid stalling 
phenomena that we encountered in computing the integrals of f% and f$. 

We point out that in the case of the function f 2 , where there exist a singularity in 
(0.4, 0.4, 0.4), triplequad uses a domain decomposition approach, implemented in the method 
quadgk (Gauss-Kronrod cubature rules) avoiding the singularity. Actually, due to the geom- 
etry of the cylinder, we could compute the exact values for the integrals by using separation 
of variables, that is instead of a call to the Matlab built-in function triplequad we used the 
product of the built-in functions dblquad and quadl. This allowed a considerably reduction 
of the computational time, as displayed in Table [6] for degrees n < 20. 

Concerning function f<±, the least-square errors seem not those that one can expected. 
In [3 Fig. 3.2] the authors already computed the relative hyperinterpolation errors for the 
Runge function w.r.t. the number of function evaluations. From that figure, correspondingly 
to polynomial degree n = 30, that requires (n + l)(n + 2)(n + 3)/6 = 5456 function evaluations, 
the hyperinterpolation error is about 10 _1 . Hence, what we see in Figure Mis consistent with 



those results and, as expected, formula (4.9) is an overestimate of the least-square error. 



n 


triplequad 


dblquad x quadl 


5 


1 min. 9 sec. 


4.7 sec. 


10 


16 min. 3 sec. 


1 min. 4 sec. 


15 


1 h 15 min. 10 sec. 


5 min. 33 sec. 


20 


3 h 47 min. 12 sec. 


16 min. 30 sec. 



Table 6: Computational time: triplequad vs doublequadx quadl 
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Figure 6: Interpolation error on AFP (blue — *), cubature errors on AFP (black — h) Interpo- 
lation error on DLP (magenta — o), cubature errors on DLP (green - ) and least-squares 
errors (red — o). The points are extracted from the WAM1. Wade basis. In abscissa the 
polynomial degree. 
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Figure 7: Interpolation error on AFP (blue — *), cubature errors on AFP (black — h) Interpo- 
lation error on DLP (magenta — o), cubature errors on DLP (green - ) and least-squares 
errors (red — o). The points are extracted from the WAM2. Wade basis. In abscissa the 
polynomial degree. 
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